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Long-term stability studies of nonlinear Hamiltonian systems require symplectic integra- 
tion algorithms which are both fast and accurate. In this paper, we study a symplectic 
integration method wherein the symplectic map representing the Hamiltonian system 
is refactorized using polynomial symplectic maps. This method is analyzed in detail 
for the three degree of freedom case. We obtain explicit formulas for the action of the 
constituent polynomial maps on phase space variables. 

Keywords: Symplectic integration; polynomial maps; Lie perturbation theory 

1. Introduction 

Numerical integration algorithms are essential to study the long term single parti- 
cle stability of nonlinear, nonintegrable Hamiltonian systems. However, standard 
numerical integration algorithms can not be used since they are not symplectic ^. 
This violation of the symplectic condition can lead to spurious chaotic or dissipative 
behavior. Numerical integration algorithms which satisfy the symplectic condition 
are called symplectic integration algorithms ^. Several symplectic integration algo- 
rithms have been proposed in the literature 2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21^ 

Some of these directly use the Hamiltonian whereas others use the symplectic map 
22,23 representing the nonlinear Hamiltonian system. For complicated systems like 
the Large Hadron Collider which has thousands of elements, using individual Hamil- 
tonians for each element can drastically slow down the integration process. One the 
other hand, the map based approach is very fast in such cases 24,25 _ 

One class of the map-based methods uses jolt factorization 6,ii,i7,i9^ g.^.^. -^jigj-e 
are still unanswered questions on how to best choose the underlying group and 
elements in the group Further, some of these methods i^'i'^'-i^ can be quite 
difficult to generalize to higher dimensions. Another class of methods uses solvable 
maps ^'^^'^^ or monomial maps Even though they are fairly straightforward to 
generalize to higher dimensions, they tend to introduce spurious poles and branch 
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points not present in the original map 

We investigate a new symplectic integration method where the symplectic map is 
refactorized using "polynomial maps" (maps whose action on phase space variables 
gives rise to polynomials). This does not introduce spurious poles and branch 
points. Moreover, it is easy to generalize to higher dimensions. Further, since it is 
map-based, it is also very fast. In this paper, we describe in detail the theoretical 
underpinnings of the polynomial map factorization of symplectic maps. We also 
apply it to Hamiltonian systems. 

2. Preliminaries 

We restrict ourselves to three degrees of freedom nonlinear Hamiltonian system. 
The effect of a Hamiltonian system on a particle can be formally expressed as the 
action of a symplectic map ^A that takes the particle from its initial state to 
its final state z^'" 2^,23 

= M 2'". (1) 
Here z represents the collection of six phase space variables: 

2 = (9l,92,93,Pl,P2,P3)- (2) 

Using the Dragt-Finn factorization theorem^^'^^, the symplectic map A4 can be 
factorized as shown below: 

M = Me'-f^'- e'-^"'- ...e'-^"'... . (3) 

Here fn{z) denotes a homogeneous polynomial (in z) of degree n uniquely deter- 
mined by the factorization theorem. The Lie transformation e''^^^^' is given by 

where 



Tl=0 



■.f{z):g{z) = [f{z),g{z)]. (5) 

Here [f{z),g{z)\ denotes the usual Poisson bracket of the functions f{z) and g{z). 
Further M gives the linear part of the map and hence has an equivalent represen- 
tation in terms of the Jacobian matrix M of the map M. ^^: 

Mzi = M.jZj = {Mz\. (6) 

The infinite product of Lie transformations exp(:/„:) (n = 3,4,...) in Eq. (3) 
represents the nonlinear part of M. 

As an application, let us consider a charged particle particle storage ring which 
typically comprises thousands of elements (drifts, quadnipoles, sextupoles etc.) Us- 
ing the above procedure, one can represent each element in the storage ring by a 
symplectic map. By concatenating these maps together using group-theoretical 
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methods we obtain the so-called 'one-turn' map representing the entire storage 
ring. The one-turn map gives the final state z^^^ of a particle after one turn around 
the ring as a function of its initial state z^^^: 

z^'^=Mz^°l (7) 

To obtain the state of a particle after n turns, one has to merely iterate the above 
mapping A'' times i.e. 

Since A4 is explicitly symplectic, this gives a symplectic integration algorithm. 
Further, since the entire ring can be represented by a single (or at most a few) 
symplectic map(s), numerical integration of particle trajectories using symplectic 
maps is very fast. 

To obtain a practical symplectic integration algorithm, we follow the perturba- 
tive approach and truncate A4 after a finite number of Lie transformations: 

M !v Me'-^^'- e-f"'- ...e^^'. (9) 

The symplectic map is said to be truncated at order P. This map is still symplectic. 

However, each exponential e'-''"' in M. still contains an infinite number of terms in 
its Taylor series expansion. We get around the above problem by refactorizing 
M. in terms of simpler symplectic maps which can be evaluated exactly without 
truncation. We use 'polynomial maps' which give rise to polynomials when acting 
on the phase space variables. This avoids the problem of spurious poles and branch 
points present in generating function methods solvable map ^^'^^ and monomial 
map refactorizations. 

3. Symplectic Polynomial Maps 

In this section we study symplectic polynomial maps in some detail. We start 
by describing the difference between monomial maps and polynomial maps with 
respect to presence of poles and branch points. This difference can be illustrated 
using the following examples. Consider the monomial symplectic map exp(: (^Tp\ :). 
Its action on gi, pi in a two dimensional phase space is given as follows: 

q'i = exp(: g^pi :)q'i = ; p'l = exp(: g'^pi :)pi = -h gi)^. (10) 

i + (Zi 

This map has a pole at gi = — 1. 

On the other hand, consider the symplectic map exp(: aiql + a2Pi :) where ai, 
a-i are real constants. We determine its action on phase space variables as follows. 
Note that the symplectic map is of the form cxp(: h{z) :) where /i(z) is a fimction 
which depends only on the phase space variables z and is independent of time i. If 
we take h{z) to be the Hamiltonian function, then solving the Hamilton's equations 
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of motion for this Hamiltonian from time t = to time t = is equivalent to the 
following symplectic map action 

z{t = t^) = exp[-(i/ - f ) : h{z) :]z{t = f). (11) 

Equivalently, obtaining the action of the symplectic map exp[— (t-'' — f) : h{z) :] on 
the phase space variables is the same as solving the Hamilton's equations of motion 
with h{z) as the Hamiltonian from time V" to . Setting = and = —Iwe have 
the following equivalence: Obtaining the action of the symplectic map exp(: h{z :) 
on phase space variables is equivalent to solving the Hamilton's equations of motion 
using h{z) as the Hamiltonian from time i = to time t = —1. In this case, z(0) 
will correspond to the initial values of the phase space variables and z{—l) to the 
final values obtained after the action of the map exp(: h{z :). 

Returning to our symplectic map, we obtain its action by first solving the Hamil- 
ton's equations of motion from time t = Q to t = —1 using the argument of the Lie 
transformation, h = a\q\ + a2Pi, as the Hamiltonian. The Hamilton's equations of 
motion are given by: 

dq\ _ dh 
dt dpi ' 

dPi _ _9h 

dt ~ dqi ^ ' 

Solving these simple equations, we obtain: 

91 (i) = 9i(0)+a2t; 

p-^(t) = pi(0)-aiait^-3aia2gi(0)i^-3aigi(0)^t. (13) 

where 91(0) and pi(0) denote the values of q\ and pi at time t = Q. To obtain the 
action of the map exp(: aigf + a^Pi :) on the phase space variables, we set t = —1 
in the above equations and denote q\{—l), p\{—l) by (/f™, p{™ and 91(0), pi(0) by 
9i"> pT respectively. Thus we get 

gf™=gr-a2, pf" = + aia^ - Saiazgi" + 3ai(gr)'- (14) 

Using Eq. (4), we can easily verify that the above result is indeed correct. We note 
that the final values of the phase space variables are polynomial functions of the 
initial variables and therefore involve no poles or branch points. This is an example 
of a polynomial map. 

We now determine the classes of symplectic maps which are also polynomial 
maps. We obtain the following simple principles which are equally applicable in 
higher dimensions. 

1. All polynomials of the form h{z) where both a phase space variable and its 
canonically conjugate variable do not occTir simultaneously give rise to 
symplectic polynomial maps via exp(: h{z) :). We will call such h{z)'s as 
polynomials of the first type. 
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2. If a canonically conjugate pair qi,pi is present in the polynomial h{z) and it 
appears either in the form [a{z)qi + g{pi,z)]"^ or [a{z)pi + g{qi,z)]"^ (where 
m = 1,2 . . ., z = {qj,Pk} with j ^ k ^ i and a, g are polynomials in the 
indicated variables), then this polynomial h{z) again gives rise to a symplectic 
polynomial map via exp(: h{z) :). If a product/sum of such factors appears 
in h{z), each term in the product/sum is a function of different canonically 
conjugate pairs. We will call h{zys of the form described above as polynomials 
of the second type. 

We can prove the above results as follows. Let z denote a collection of phase 
space variables {qj,Pk} with j ^ k. Thus polynomials of the first type are of the 
form h(z). The polynomial map is then given by exp(: h{z) :). As described earlier, 
its action on the phase space variables is given by solving the Hamilton's equations 
of motion from time t = to t = —I using h{z) as the Hamiltonian. From classical 
mechanics each of the variables in the collection i is a cyclic variable and is 
therefore conserved by the Hamiltonian. Thus z{t) = z{0). Consequently 

f/*" = exp(: h{z) :}z"'' = F". (15) 

Next we consider the action of this map on the variable Zj which is canonically con- 
jugate to Zj. Solving Hamilton's equations of motion with h{z) as the Hamiltonian 
from time t = to t = —1 we get 

£,-(-!) = f,(0)-(-l)'-^(z,(0)), (16) 

where r is zero if Zi is a coordinate variable and 1 otherwise. As before, the action 
of exp(: h{z) :) on Zj is obtained by setting Zj{0) as fj*" and ij (— 1) as f^-^"': 

^/" = i'/"-(-ir^(^r)- (17) 

Since h{z) is a polynomial in z, the right hand side of the above equation is a 

polynomial in z and z. From Eqs. (15) and (17) we conclude that exp(: h{z) :) is 
a polynomial map. Further, since all Lie transformations are symplectic maps 
exp(: h{z) :) is a symplectic polynomial map. 

Next we consider polynomials h{z) of the second type described in item 2) above. 
Let h{z) — [a{z)qi + g{pi, z)]™ where qi,Pi is a canonically conjugate pair. Further 
m is a positive integer, z = {qj,Pk} with j ^ k ^ i and a, g are polynomials in the 
indicated variables. We will show that exp(: h{z) :) is a polynomial map. The proof 
for the case h'{z) = [a{z)pi + g{qi, z)]™ is similar. A concrete example of the type 
of h{z) that we are considering is given by h{z) = {a\P2 + 0:2^3 + ctsp^ + 04^2^3 + 
oibPz + Q!693P2)™ where a^'s are real constants. 

As before we first solve the Hamilton's equations of motion with h{z) as the 
Hamiltonian. Since each of the variables in the collection denoted by z is cyclic, we 
have z{t) = 2(0) and 

z^'"- = exp(: h{z) :)^™ = (18) 
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Solving the Hamilton's equation of motion for pi we get 
dpi dh{z) 



dt dqi 



m[a(z(0))ft(0) + g{pi{0), z{0))r~'a{m), (19) 



where we have used the fact that [a{z)qi + g{pi, z)] is a conserved quantity under 
this Hamiltonian flow i.e. [a{z)qi+g{pi,z)] = [a{z{0))qi{0)+g{pi{0),z{0))]. Solving 
this simple equation we obtain 

Pi{t) = Pi{0) - m[a{z{0))qi{0) + g{pi{0), z{0))r-'a{z{0))t. (20) 

Note that Pi(t) is a polynomial in both z and t. Setting t — ~l and denoting 
z(0),z(— 1) by z^^-jZ^^" respectively, we see that pf*" = exp(: h{z) •.)pl^ is a 
polynomial in z'": 

" = Pf + m[a(r-)C) + 9{pT, z^nr-'a{z'^). (21) 
The equation of motion for qt gives us 

f = ^ = -m[amMo)+Mo),mr-'^^. (22) 

Since Pi(t) is a polynomial in z and t, the right hand side of the above equation 
is also a polynomial in z and t. Therefore, the above differential equation can be 
easily integrated to give qi{t) which is guaranteed to be a polynomial in z and t. 
Setting t = —1 and denoting 2;(0),2;(— 1) by ^™^^/»" respectively, we obtain the 
result that g/'" = cxp(: h{z) :)(jrj" is a polynomial in z'". Finally, the equation of 
motion for zj, the variable canonically conjugate to zj, is given by 



dt = ^"^^'^ = {-lYm[a{z{0)MO) + 9{Pi{0), m)Y 

d[a{z)qi+ g{p^,z)] 



dzj 



-iz = m), (23) 



where r is zero if zj is a coordinate variable and 1 otherwise. The right hand side is 
a polynomial function of qi,Pi both of which in turn are polynomials in t. Hence the 
equation can be integrated giving Zj{t) as a polynomial in z and t. Setting t = —1 
and denoting z{0),z{—l) by respectively, we find that Zj-^™ = exp(: h{z) : 

)zj" is a polynomial in z'". Thus we have proved that exp(: h{z) :) is a (symplectic) 
polynomial map. 

To conclude, we consider the case where a sum/product of factors of the form 
[a{z)qi+g{pi, z)]™ or [a{z)pi+g(qi, z)]"' appear in h{z). If a phase space variable ap- 
pears in one factor of the sum/product, by assumption, neither this variable nor its 
canonically conjugate variable appears in the remaining factors of the s\im/product. 
Therefore each term in the sum/product is independently conserved by the Hamil- 
tonian flow generated by h{z) and acts only on the phase space variables appearing 
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in that term. Consequently each term can be considered separately and since each 
term is of the form [a{z)qi + g{pi, z)]™ or [a{z)pi + g{qi, z)]™, the argument given 
above can be immediately applied proving that exp(: h{z) :) is a symplectic poly- 
nomial map even in this case. This completes the proof of the claims made in items 
1) and 2). We conjecture that all symplectic polynomial maps have one of the two 
forms enumerated above. 

4. Symplectic Integration using Polynomial Maps 

In this section, we return to the problem of symplectic integration. We restrict 
ourselves to symplectic maps in a six dimensional phase space truncated at order 
4. The results obtained below can be generalized to both higher orders and higher 
dimensions using symbolic manipulation programs. The Dragt-Finn factorization 
of the symplectic map is given by: 

M = Me-f^'e-^'''-, (24) 

where 

/s = a2sql + o.2miPi H H asap!, 

fi = as4qt + asaqfpi H h 0209^3. (25) 

Here the coefficients a28, ■ ■ ■ , 0209 can be explicitly computed given a Hamiltonian 

system^^ and arc therefore known to us. The numbering of these monomial coeffi- 
cients follows the standard Giorgilli scheme . The above map captures the leading 
order nonlinearities of the system. Since the action of the linear part M on phase 
space variables is well known [cf. Eq. (6)] and is already a polynomial action, we 
only refactorizc the nonlinear part of the map using iV polynomial maps This 
is done as follows: 

M!=^r = Me''''-e-^^'----e-^'''-, (26) 

where e'^'^^'s are symplectic polynomial maps and the numeral appearing in the 

subscript indexes the polynomial maps. The polynomial maps are determined by 
requiring that V agree with M up to order 4. That is, when the N polynomial 
maps are combined, the resulting symplectic map should have all the monomials 
present in /a and /4 with the correct coefficients up to order 4. 

The basic idea in obtaining the required refactorization is to group the monomial 
terms present in /a and [cf. Eq. (25)] such that the Lie transformation corre- 
sponding to each grouping gives a polynomial map. It is obviously easy to handle 
monomials where both members of the canonically conjugate pair are not present 
simultaneously. They can be grouped into different polynomials (for example, mono- 
mials involving only the coordinate variables g^'s in one group and those involving 
only the momentum variables p^'s in another group etc.) so that each one of these 
is a polynomial of the first type. Since a product of two Lie transformations (e''^^' 
and e■^*^) is being refactorized as a product of many simpler polynomial symplec- 
tic maps, the coefiicients multiplying the monomials in each individual polynomial 
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map will be in general different from the coefficient multiplying the corresponding 
monomial in Eq. (25). The relation between these coefficients is easily obtained 
using the CBH theorem Monomials s{z) where both members of the canonically 
conjugate pair are present simultaneously (like for example, s{z) = QiPi) are more 
difficult to handle since the corresponding Lie transformations e'®*^^)' typically give 
rise to poles and branch points which we wish to avoid. But by using a product of 
two polynomial maps of type 2 with carefully chosen coefficients, these can also be 
generated. 

Using the above procedure, it turns out that we require 23 polynomial maps for 
refactorization; 

M^V = Me-''''- e-''^'- ■ ■ ■ e^'^^^^ , (27) 
The hiS are given as follows: 

hi = qf 628 + qf 92 &30 + 9i 93 &32 + 9i 92 ^39 + 9i 92 93 hi + 9i 9.3 he + 

ql hi + ql 93 h& + 92 93 hi + 93 ^80 + 9i ^84 + 9i 92 ha + 9i 93 ^88 + 
9i 92 hb + 9i 92 93 hi + 9i 93 ^102 + 9i 92 ^120 + 9i 92 93 6122 + 

91 92 93 ^127 + 91 93 ^136 + 92 ^175 + 92 93 '>177 + 92 93 ^182 + 

92 93 ^191 +93^205, 

h2 = [(&29 + hi) + 92 {hi + '>i06) + P2 {h2 + hor) + 93 {hs + hos) + 

P3 (&94 + &IO9)] (Pl + 91)^, 

ha = [(-&29 + ^34) + 92 (-691 + hoe) + P2 (-692 + ^107) + 

93 (-^93 + fel08) + P3 (-^94 + &I09)] {-Pl + 9l)^, 

hi = [{hb + hs) + 91 (^121 + ^124) + Pi ihbG + ^159) + 

93 {hso + hse) + P3 {hsi + hs?)] {P2 + 92)^, 
hs = [{-hb + hs) + 91 {-h2i + &124) + Pi (-&156 + &159) + 

93 (-&180 + hse) +P3 {-hsi + hs?)] {-P2 + 92)^, 
he = [{hi + h2) + 91 {hsr + hss) + Pi (&172 + &173) + 92 {h92 + hm) + 

P2 {ho2 + ho3)] {P3 + 93)^, 
hr = [{-hi +^82) + 91 (-^137 + ^138) + Pi {-h72 + hva) + 

92 (-^192 + ^193) +P2 (-&202 + &203)] {-P3 + 93)^, 



hg 


= {Pi + 


91) 


(92 hb - 


1- 93 ^37 H 


- ql hio - 


f 92 93 ^'112 


+ P3 92 ^113 + 93 hn) 


hg 


= {Pi + 


91) 


{P2 he - 


f Pz ha - 


\-pl hii 


+ P2 93 ^115 


+ P2P3 hie +pI hi9 


hio 


= {P2 + 


92) 


(91 ha - 


'r 93 &69 H 


- ql he + 


91 93 &125 ^ 


-Psqi &126 + 9.3 ^188) 


hii 


= {P2 + 


92) 


{pi hh - 


t- P3 ho - 


\- pi hie ■ 


+ Pi 93 ^160 


+ P1P3 hei +pI hgo 


hi2 


= (P3 + 


93) 


(91 ^47 - 


\- 92 ^^72 H 


- 9l ^103 " 


f 91 92 ^128 


+ P2 91 ^134 + 92 ^183) 


hi3 


= {P3 + 


93) 


{pi h2 - 


1- P2 he. - 


f p\ hb3 


+ Pi 92 ^163 


+ Pi P2 h69 + pI h99 


hii 


= P2 Ql hi + Ps 9i h3 + pI 9i h3 + P2 P3 9i hb + P3 9i ^48 + P2 ql h? + 
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P3 qf b89 + pI 1i ^99 + P2 P3 ll &101 + pI 5i '>104 + pI 91 &130 + 

pI P3 qi bi32 + P2 pI qi &135 + pI qi ^139 , 
hi5 = p\ q2 &50 + pI qs &52 + Pi ql hi + Pi 92 93 he + Pi ql hi + p\ 92 bm + 
p\ qs &143 + pI ql &145 + p\ q2 <?3 &147 + pi ql ^'152 + pi ql &155 + 
pi ql <?3 ^157 + pi q2 ql bie2 + pi ql bm, 

hl6 = Pi P3 q2 b57 + P3 ql b67 + pj q2 ^73 + Pi P3 «2 ^148 + Pi P3 ql ^158 + 
Pi pI 92 &164 + P3 ql bl78 + pI ql &184 + pI ?2 6l94, 

hn = P2 qi qs ba + pi 93 675 + p2 ql 677 + P2 ql qs &100 + pI qi 93 &131 + 

P2 91 93 ^133 + pI 93 ^196 + pI ql bl98 + P2 ql 62OI, 
his = Pi P2 93 &59 + J32 93 ?>150+ Pi P2 93^166+ Pi P2 93 ^168, 

hl9 = P3 9l 92 ^42 +P3'?1 92 &98 +P3 91 92 ^123 +P391 92&129, 
h20 = P?&49 +P?P2&51 +P?P3&53 +P1P2&58 +P1P2P3&6O +P1P3&63 + 
pI 674 + pI Pa ^76 + P2 pI &79 + pI bs3 + p\ &140 + Pi P2 ^142 + 
P\ P3 buA + pI pI fcl49 + pI P2 P3 &151 + P? pI &154 + Pi pf &165 + 
Pi P2 P3 &167 + Pi P2 P3 &170 + Pi P3 ^174 + P2 ^195 + P2 P3 &197 + 
P2 Ps ^200 + P2 Ps &204 + P3 &209, 
h21 = (Pl +91 +Pl&105)^ + (P2 + 92 +P2^'185)^ + (P3 + 93 +P3&208)^, 

h22 = {-Pi - qi + qlbsbf + {-P2 - q2 + qlbnof + 

(-P3 - 93 + 93 ^206)^, 

4 2 2 2 2 

h23 = (pi+9i) &90 + (Pi+9i) (P2 + 92) &iii + (pi+9i) (P3 + 93) &118 

+ (P2 + 92)^ &179 + (P2 + 92)^ (P3 + 93)^ &189 + (P3 + 93)^ &207- 

Here &i's are at present unknown coefficients. As mentioned above, by forcing the 
refactorized form V to equal the original map M. up to order 4 and using the 

CBH theorem^^, we can easily compute these unknown coefficients in terms of the 
known a^'s. These expressions are available from the author as part of a FORTRAN 
program implementing the above algorithm. 

The explicit actions of the polynomial maps on phase space variables can be 
obtained and they are given below. This completely determines the refactorized 
map V. Each exp(: hi :) is a polynomial map which can be evaluated exactly and 
is explicitly symplectic. Thus by using V instead of M. in Eq. (8), we obtain an 
explicitly symplectic integration algorithm. Further, it is fast to evaluate and does 
not introduce spurious poles and branch points. The above factorization is not 
unique. However, the principles outlined earlier impose restrictions on the possible 
forms and this cases considerably the task of refactorization. Moreover, we require 
the coefficients bi to be polynomials in the known coefficients . Otherwise this can 
lead to divergences when a^'s take on certain special values. Finally, we minimize 
the number of polynomial maps in the refactorized form. Our studies show that 
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different polynomial map refactorizations obeying the above restrictions do not lead 
to any significant differences in their behavior. 

We now derive the explicit actions of the polynomial maps on phase space vari- 
ables. First consider exp(: hi :). We obtain its action on the phase space variables 
by following the procedure outlined in the paragraph before Eq. (12). Wc notice 
that hi = hi{qi,q2, qs) depends only on the coordinate variables which are therefore 
cyclic variables. Hence we immediately obtain: qi{t) = qi{0) {i = 1,2,3). Solving 
the Hamilton's equations of motion for pi with hi as the Hamiltonian from t = 
to t = —1 we get: 

Bh 

Pi(-l)=Pi(0) + ^((?i(0), 92(0), 93(0)), i = 1,2,3. (28) 

Denoting Zi{0) and Zi{—1) by zl" and z{^" respectively, we finally obtain the action 
of exp(: hi :): 

gf" = 9r; pf"=pr + ^(3r,9r,9r), * = 1,2,3. (29) 

Next consider the action of exp(: h2 ■)■ The 2 factors in the product are inde- 
pendently conserved under the Hamiltonian flow. Consequently, we have 

A2 = [(629 -I- 634) + 92 (&91 + bioe) + P2 (692 + &107) + 

93 (&93 + ^lOs) + P3 (^94 + &IO9)] 

= [(&29 + bu) + 92(0) (691 + bloe) + P2{0) (692 + &107) + 

93(0) (693 + bios) +P3{0) (&94 + biog)] , (30) 

B2 = (9i+Pi) = (9i(0)+pi(0)). 

Solving Hamilton's equation of motion with /12 as the Hamiltonian from t = to 
t = —1 we get 

qi{-l) = gi(0)-3^2B|; pi{-l) = pi{0) + 3A2BI 

92(-l) = 92(0) -53(692 + 6107); P2i-1) = P2{0) + Bl {hi + bloe) ; i^l) 

93(-l) = 93(0) -53(694 + 6109); P3(-l) =P3(0) + 53(693 + 6108). 

Denoting 2,(0) and Zi{—1) by z}"" and 0/'" respectively, we obtain the action of 
exp(: h2 ■■)■■ 

q('^ = ql--3A2Bl; pf " = + 3^25^; 

= 9^-51(692 + 6107); pr=pr + 5|(69i+6io6); (32) 
qi'^ = <zr - 52^ (694 + 6109) ; pt =pT + Blib,, + bios), 

where A2, B2 are now functions of z*". The actions of cxp(: hi :). i = 3, 4, . . . , 7 
on the phase space variables are obtained in a similar fashion and these actions are 
listed in the Appendix. 
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We now consider the action of exp(: hg :). From the Hamilton's equations of 
motion, we have the foUowing conserved quantities: 

^8 = (92 &35 + 93 ''ST + 92 ^110 + 92 93 ^'112 + PS 92 &113 + 93 ^ll?) 

= (92(0)635 + 93(0)637 + 92(0)6110 + 92(0)93(0)6112+ 

^3(0)92(0)6113 + 93(0)6117), (33) 
Bs = (9i+Pi) = (gi(0)+pi(0)). 

Solving the equations of motion for 9j's with /ig as the Hamiltonian from t = to t 
we get: 

qi{t) = qiiO) + 2AsBst; 92 (i) = 92(0); 

93(t) = 93(0)+B|6ii392(0)i. (34) 

For the momentum variables we get the following differential equations: 
dpi{t) 



dt 
dp2{t) 

dt 
dpajt) 

dt 



-2AsBs, 

-Bl[bs5 + 26iio92(t) + 6ii293(i) + bnaPait)], (35) 

1292(0 + 2611793 (t)]. 



The first equation can be trivially solved to obtain pi{t). After substituting for 
92 (t) and qsit) which are known, we can next solve the last equation for p^it). 
Substituting this in the second equation, we finally get P2{t)- Setting t = —1 and 
denoting -2^(0), Zi{—1) by zf^^ respectively, we obtain the action of exp(: hs ■)■ 

gf- = qY^-2AsBs; p(''' = pT + 2AsBs; 

9I" = qT; P^'"=Pr + S|[635 + 26iio9r + &ii29r + 6irf+ 

B|6n3 (637/2 + 6ii39f - -B8'''ii36ii79r/3)] ; (36) 
9I'" = qt-BlWmT\ 

pf'^ = pf + Bl [637 + 6ii29r + 26ii79f - S|&ii36ii79r] , 

where Ag, B^ are now functions of 2;'". The actions of exp(: hi :), i = 9, 10, . . ., 13 
on the phase space variables are obtained in a similar fashion and these actions are 
listed in the Appendix. 

Next consider the action of exp(: hu :). We notice that /114 = hu{q\,p2,Pa) is 
independent of pi, 92 and 93. Hence 91, P2, Ps ai'C' cyclic variables and are conserved 
under the action of the Hamiltonian. Solving the Hamilton's equations of motion 
for pi, 92 and 93 with hu as the Hamiltonian from < = 0to< = — Iwe get: 

pi(-l) = pi(0) + ^(9i(0),p2(0),P3(0)), 
C91 
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<Z2(-1) = 92(0) - ^(gi(0),p2(0),i>3(0)), (37) 
qsi-l) = <?3(0) - ^(gi(0),p2(0),P3(0)). 

<JP3 

Denoting Zi{0) and Zi{—1) by z}" and 0/'" respectively, we finally obtain the action 
of exp(: hi4 :): 

oqi 

rih 

fin _ in '-"nA f jn in in\ ^fin _ in 
OP2 

^fin ^ pf ^ PT ■ 

The actions of cxp(: hi :), i ~ 15, 16, ... , 20 on the phase space variables are ob- 
tained in a similar fashion and these actions are listed in the Appendix. 

We now consider the action of exp(: /121 :)• From the Hamilton's equations of 
motion, we have the following conserved quantities: 

A21 = (pi + gi+p?6io5) = (pi(0) + 91 (0)+pf (0)6105), 

B21 = {P2 + 92 +P2 &185) = (P2(0) + 92(0) +pI{0) 6185) (39) 

C21 = {P3 + 93 + P3 ^208) = (P3(0) + 93(0) + pI{0) 62O8) • 

Solving the equations of motion for pi's with /121 as the Hamiltonian from t = to 
t we get: 

pi(t) = pi(0)-3Aiit; p2(i) =P2(0) -3B|ii; 

Piit) = p3(0)-3C|ii. (40) 



The equations of motion for gj's are given as: 

= 3Aii(l + 26io5Pi(t)) 



<^9i(*) _ o a2 



dt 

= 3B2'i(l + 26i85P2(t)), (41) 



'^92(0 ^„2 



dt 

^ : 352^(1 + 2fe208P3(i))- 



rf93(i) _ od2 



dt 

Substituting the expressions for pi (t) , p2 (^) , Pz (t) obtained earlier, the above equa- 
tions can be easily solved. Setting t = —1 and denoting Zi{0), Zi{—1) by zj", z/*" 
respectively, we obtain the action of exp(: /121 :): 

gf" = ql"-3AUl + 2bio5PT + -'iAl,bw5); p{"' =pT + 'SAI,; 
qi'^ = qi"-3Bl,il + 2bis,pt + Wlihs5); p^ ^ pT + 'iBi,; 
qt = 9r-3C|i(l + 26208Pr + 3C|i6208); pi'"" = pf + ^C^^, 
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where ^21,-621,(721 are now functions of 2;™. The action of cxp(: /122 on the 
phase space variables is obtained in a similar fashion and is listed in the Appendix. 

Finally, we consider the action of exp(: /i23 :)■ From the Hamilton's equations 
of motion, we have the following conserved quantities: 



(42) 



A23 


= (pi + qi) 


= (pi(0) + gi(0)). 


B23 


= {p2 + q2) 


= b2(0)+(72(0)), 


C23 


= iP3 + 93) 


= (P3(0) + (Z3(0)). 



From the equations of motion with /123 as the Hamiltonian we get: 



— — — 4O90A23 + Z0iiiA23±i23 + ^Oll8^23(^; 

dq2{t) _ dp2{t) _ jo3 .r,, 42 R ,9, R ^2 

— 40i79i323 + ^0111^23 -°23 + ^Oi89-D23L^23 



dt dt 

^ = 46207CI3 + 26118^^3(^23 + 26189B23C23. (43) 



dq3{t) _ dpsit) _ ^3 , 01, a2 ry , oj. d2 



dt dt 

Solving these equations from t = to t = —1 and denoting Zi{0), Zi{—1) by 
2:/*" respectively, we obtain the action of exp(: /123 :): 

qfin ^ ^ra _ [4,,gQ^3^ + 26iii A23SI3 + 26118A23CI3] i 
pfin ^ ^ [4j,g^^3^ ^ 26111 A23BI3 + 26ll8^23C|3] \ 

g|™ = ^2" ~ [46179SI3 + 26111^23-^23 + 26189B23CI3]; 

Pt = + [46179^13 + 26111^13^23 + 26189B23CI3]; (44) 

13" = 'it - [46207CI3 + 26118^23^23 + 26189-BI3C23]; 

p/™ = Pf + [46207CI3 + 26ii8^i3C23 + 26i89S2'3C23], 

where A23 , -B23 , C23 arc now functions of . 

Substituting in Eq. (27) the explicit formulas for the actions of the polynomial 
maps listed above and in the Appendix, we can evaluate the action of V without 
violating the symplectic condition. Using this explicitly symplectic map in Eq. (8), 
we have the desired symplectic integration algorithm. 



5. Applications 

We have applied the method to a large particle storage ring for storing charged 

particles. This storage ring consists of 5109 individual elements (where these ele- 
ments could be drifts, bending magnets, quadrupoles or sextupoles). If one tries to 
numerically integrate the trajectory of a charged particle through this ring using a 
conventional integration algorithm, one has to go through the ring element by ele- 
ment where each clement is described by its own Hamiltonian. This is cumbersome 
and slow and further, does not respect the Hamiltonian nature of the system. On 
the other hand, a map based approach where one represents the entire storage ring 
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in terms of a single map is much faster ^^i^^. When this is combined with our poly- 
nomial map refactorization, one obtains a symplectic integration algorithm which 
is both fast and accurate and is ideally suited for such complex real life systems. 
The qi — p\ phase plot for one million turns around the ring using our polynomial 

map method is given in Figure 1. In this case, qi and pi represent the deviations 
from the closed orbit coordinate and momentum respectively. From theoretical con- 
siderations, we expect the so-called betatron oscillations in these variables. This 
manifests itself as ellipses in the phase space plot of qx and pi variables. In Figure 
1, we observe the expected betatron oscillations. We also see the thickening of the 
ellipses caused by nonlinearities present in the sextupoles. 

0.0004 I 1 1 1 1 1 1 1 

0.0003 - 

0.0002 - 

0.0001 - 

D. - 

-0.0001 - 

-0.0002 - 

-0.0003 - 

-0.0004 I 1 1 1 1 1 1 1 

-0.0004 -0.0003 -0.0002 -0.0001 0.0001 0.0002 0.0003 0.0004 

Fig. 1. This figure shows the q\ — pi phase space plot for one miUion turns around a storage ring 
using the polynomial map method (only every 1000th point is plotted). 

6. Conclusions 

To conclude, we described in detail a new symplectic integration algorithm based 
on polynomial map refactorization. We enumerated the types of symplectic maps 
which give rise to polynomial actions on phase space variables. For a six dimensional 
phase space, we obtained the refactorization of a given symplectic map in terms of 23 
polynomial maps. The explicit actions of these polynomial maps were derived. This 
polynomial map method can be used to study long term stability of complicated 
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nonlinear Hamiltonian systems. 
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Appendix A 

The actions of polynomial maps exp(: hi :) which were not listed in the main text 

are given in this Appendix. 

The action of cxp(: /13 :) is given as follows: 

q^'^ = qT-Bl{-b,2 + bio7); p|*"=pr + Sf (-&91 + &106); (1) 

= g»-5|(-694 + 6l09); pf"=Pr + ^l(-^'93 + 6l08), 

where 

^3 = [(-^29 + 634) + qT {~bgi + feloe) + pT (-692 + &107) + 

93" (-^93 + bios) + pT {-b94 + 6109)] , (2) 

B3 = (qr-pf)- 

The action of exp(: /14 :) is given as follows: 

= C - ^I(&i56 + &159) ; pf" = K" + ^1(^121 + &124) ; 

gl'" = 9^-3^454'; pt=PT + 3A,Bl; (3) 
= qT-Blibisi + bis7); pi'"" ^pT + Bl{biso + bis6), 

where 

A4 = [(665 + bes) + qT (bi2i + 6124) + pT (^156 + ^'159) + 

qT (biso + 6186) + pT (bisi + bisr)] , (4) 
B4 = (qT+PT)- 
The action of exp(: /15 :) is given as follows: 

gf" = gi" - S| (-^156 + &159) ; pr"=pf + SIM121 + 6124); 

qt = cfr + ^A,Bl- pi'"^pl- + 3A,Bl, (5) 

g/- = q^-Bli-hsi + bisr); pt = pT + BI {-biso + bise) , 



where 



A5 = [(-665 + 668) + qi" (-6121 + 6124) + pT (-6i56 + 6159) + 

qT (-&180 + 6186) + pT (-^-ISI + 6187)] , (6) 
B5 = {qT-pf)- 
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The action of exp(: he :) is given as follows: 

gf" = gi" - B|(6i72 + 6173) ; pC" =pt + Bl{h37 + h3s); 

gf" = gr - ^6 (^202 + &203) ; P2^"=Pr + S|(&192 + bl93); 

where 

^6 = [(&8I +bs2) + qT {h37 + bias) + pf {bm + &173) + 

C (^192 + &193) (^'202 + ^203)] , 

B6 = (qT+PT)- 
The action of exp(: hr :) is given as follows: 

gf" = gf - ^7 (-&172 + &173) ; pC'' =pT + BU-bl37 + bi3s); 

qi'^ = qT-BU-b202 + b203); pt ^ pT + B^-b^m + bx^s) ; 
qt = qT + ^ArBl p^'^ = pf + 3ArBl 

where 

A7 = [(-681 + ^82) + qf {-bi37 + hi?.^) + pf (-&172 + &173) + 

qT {-bi92 + bi93)+pT (-^202 + &203)] , 
B7 = (93 +-P3 )■ 
The action of cxp(: hg :) is given as follows: 
gf" = qr-2AgBg; p{™ = pi" + 
g/'" = g™ + Bj [-536 " 26114?^ - bn5qT ~ bn6PT+ 

Blbn5{b38/2 + bn9PT + BibnsbngpT /3)] ; = pT 

g|*" = qT - Bl [638 + 6ii6Pr + 26ll9Pr + S9^115&119Pr] ; 

where 

A9 = {pTb3,+ptb3S + {pTfbru+ptqTbii5+pTptbii6 + 

ipTfbiw), 
B9 = (qT+pT)- 
The action of exp(: hio :) is given as follows: 

gf" = "= + ^fo [^40 + 269691" + &125«r + ^1" + 

Bfabi26{be9/2 + b,ssqT ' S?o&i266i88grV3)] ; 
g|'" = gr-2ylioSio; = P2" + 2^10^10; 

g|- = qt-Blobi26qT; 

pi'" = pT + Bio [b69 + 612591" + 261889^ - Si'o6i266i88gr] , 
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where 



^10 = {qr ho + qt h9 + {qTf &96 + qT qt bi25 + pf gj" &126 

HqtfhBi) (14) 
Bio = {qT+pt)- 

The action of exp(: hn :) is given as follows: 

q('" = ql" + Bf, [-655 - 26i46Pr - bieoqT " bieipt+ 

Bf,bieo{b7o/2 + bigopT + Snfoi6o6i9oPr/3)] ; pf" = pTi 
= 4"-2^nSii; pf"=pr + 2^iii?ii (15) 
= qT - Bfi [bro + bieipt + '^bi^opf + Bl^buobisopT] ; 



pf" = pt + Bf,bieoPT, 



where 



All = {pT b55 + pT bro + {pff &i46 + pT qT bi60 + pT pT &161+ 

(pr)'M' (16) 

fill = (qT+PT)- 
The action of exp(: hi2 :) is given as follows: 



where 





= qf; pt=pT + Bl^ 


[b47 + 26i039r + ^1289^ + bl34PT + 




Bl2bMb72/2 + bis3qT 


- Bl^bis4bis3qT/S)] ; 


qt 


= q™ ~ i??2^i349r; 






= pT +Bl^[b72 + bi2&qT 


+ 2bis3qT - Bl^biubis^qT] ; 


qt 


= qr-2Ai2Bi2; pt = 


--pf + 2Ai2Bi2, 


A12 


= {qTbir + qtb72 + {qT 


fbio3 + qrqtbi2s+pTqrbi34+ 




[qTfbiss) 




B12 


= {qt+pf}- 





(17) 



(18) 



The action of exp(: /113 :) is given as follows: 

gf» = q{^ + Bf, [-662 - 2bi53PT - bmsqT - bl69PT+ 

Si's^-ies (678/2 + 6i99Pr + Bf,bie3bi9opT/3)] ; pf" = pf; 
qt = qT - Bfs [b78 + bie9PT + 2^199^^ + Bl^biesbwBPT] ; 
pt = pT + Bf,bie3pT (19) 
qt = qr-2Ai3Bis; pt = pT + 2AisBis, 
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where 

Bi3 = {qT+pT)- 

The action of exp(: hi5 :) is given as follows: 

^fin ^ qY'-^{pt,qt,'ity, /r=pt, 

qt = 4"; pt=PT + ^{pT,qT,qT), 

oq2 

qt = qt; pt=pt + ^{pt,qT,qt)- 

oqs 

The action of exp(: /iig :) is given as follows: 

qt = qT; pi'''=PT + ^ipt,qT,pt), 

qt = qT-^ipT,qT,pT); pi"'=pT- 

OP3 

The action of exp(: hn :) is given as follows: 

„fin _ in. fin _ in , ^^17 / m in in\ 

q\ — «i ) P\ —Pi +^ v9i >P2 ^ 

aqi 

92 — 92 S yll )P2 )93 J) P2 — P2 J 

OP2 

Jin _ „m _ in , ^^17 / m in in\ 

93 — 93 I P3 — P3 + o„ V9l 1P2 )93 j)- 

0'93 

The action of exp(: hi% :) is given as follows: 

^/i„ ^ gin_^^^»^^»^^»). p/-=pin^ 

qt = qT-^ipT,pt,qt), pt=pf, 

OP2 

rlh 

93 — 93 7 -P3 — -Ps + IPl )P2 )93 >■ 

The action of exp(: /iig :) is given as follows: 

9^" = 9r; p^r=pT + ^{qT,qT,pt), 

qt = qT; Pt =PT + ^{qT,q2\pr), 

oq2 

„fin _ in 9hig ■ ■ ■ fin _ in 

93 ~ 93 v9i )92 )P3 /' P3 —Pa ■ 

<jpi 



Polynomial Map Symplectic Algorithm 



The action of exp(: /120 is given as follows: 

qt = qT-^{pt,pt,pt), Pt=pf, (26) 

Jin _ „m _ t^"-20 , in in in\ tin __ in 

93 — 93 \Pl 5^2 '^3 J) -P3 — P3 ■ 

The action of exp(: h22 '■) is given as follows: 



= 


qT 


+ 3^22 ; 




= pT 


- 3^^2(1 


- 268531" - 


- 3A22685); 


qt = 


qT 


+ 3-BI2; 


P^r 


= pt 


- 35^2 (1 


- 26i769r 


— 3B|2^176); 


qt = 




+ 3CI2 ; 


Pt 


= pt 


- 3C|2(1 


— 2620693" 


— 3CI262O6), 



where 



(27) 



^22 = ={-pT'qT + {qT)'bs5), 

B22 = {-pT-qT + iqTfh76) (28) 
C22 = {-PT-qT + iqTf hoe)- 
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